Is there a Reproducibility crisis?

Main classes of ODE models, ordered by increasing complexity¹

• Linear Equations
• Michaelis-Menten
• Homoeostatic, Regulatory Models

• Constant fluxes, and/or Production and decay model.

• Example: secretion of IL-2 by Th1 cells minus IL-2 decay, in Equation 1:

  • \(I_2\): cytokine concentration
  • \(T_1\): substrate, cell-type abundance
  • \(\nu_{21}\) and \(\delta_2\): parameters of the model, respectively production rate and decay rate

\[ \frac{dI_2}{dt} = \nu_{21} T_1 - \delta_2 I_2 \qquad(1)\]

• Catalytic activation: used when a the speed of the reaction is driven by a limiting resource:

• Lo et al. (2016) models reactions’ kinetics with a Michaelis-Menten process.

• Example: differentiation of \(M_0\) macrophages into \(M_1\) functional macrophages driven by \(TNF-\alpha\), see Equation 2:

  • \(f_1\): (optional), basal activation of \(M_0\)
  • \(M_0\): undifferentiated substrate
  • \(I_{\alpha}\): modifier or catalyst of the reaction (e.g., antigen, cytokine precursor)
  • \(\sigma_{M_{\alpha}}\): maximum rate signalling or activation rate.
  • \(\zeta _{M_{\alpha}}\): half-saturation constant

\[ \frac{dM_1}{dt} = \left(f_1 + \frac{\sigma_{M_{\alpha}} I_{\alpha}}{\zeta _{M_{\alpha}} + I_{\alpha}}\right) M_0 \qquad(2)\]

• Uncompetitive inhibition, formula is the inverse of the catalytic activation.

• Example: alteration of IL-\(21\) activation efficiency by the IL-4 interleukin, in Equation 3:

  • \(M\): the macrophages subset
  • \(\zeta_4\): Michaelis half-constant for inhibition activity
  • \(I_{21}\) and \(I_4\): Activator and Inhibitor concentration, respectively

\[ \frac{dT_{17}}{dt} = M \frac{1}{1 + \frac{I_4}{\zeta_4}} \qquad(3)\]

• Complex and intertwined network of positive and negative Feedback loops.
• Example: Differentiation of T-17 population is primarily driven by macrophages induction.
• However, the overall efficiency is regulated by a complex balance of cytokines activation and repression, see Figure 1:

All ODE equations adhere to the rate-law principle, see Tip 1.

2) Sensitivity analyses

ii) Latin HyperCube Sampling, step I

Main parameters:

• \(x_i^{(j)}\) is the \(i\)-th sample drawn for the \(j\)-th dimension

• \(d=9\) (number of parameters included in the sensitivity analysis)

• \(N=5000\) the number of parameter distributions simulated

• \(\mathcal{U}(x_{\min_j}^{(j)}, x_{\max_j}^{(j)})\) is the Uniform distribution, with bounds defined independently per dimension.

2) Sensitivity analyses

ii) Latin HyperCube Sampling, step II

n_samples <- 5000

parameters_sensitivity <- parameters_healthy$value
param_ranges <- tibble::tibble(parameters_name = parameters_healthy$key,
                               LB= parameters_sensitivity *0.8,
                               UB = parameters_sensitivity *1.2)


set.seed(20)
bootstrap_parameters <- lhs::randomLHS(n_samples, length(parameters_sensitivity))
for (i in seq_along(parameters_sensitivity)) {
  bootstrap_parameters[,i] <- qunif(bootstrap_parameters[,i],
                                    param_ranges$LB[i],
                                    param_ranges$UB[i])
}
apply(bootstrap_parameters,2, quantile)

1

Determine hyper-parameters of sensitivity analysis, such as the number of simulations \(N\), and lower and upper bounds within a \(\pm 20\%\) range.

2

We fix the random seed to ensure reproducibility of the sensitivity analyses.

3

The core function of the lhs package, lhs::randomLHS(), is only able to perform standard uniform sampling.

4

Apply the Inverse Transform Sampling Theorem (see Tip 2) to switch from a standard uniform distribution (\(U \sim \mathcal{U}(0, 1)\)) modelling to any bounded uniform distribution¹

5

Apply the quantile function to verify the quantiles are rather close from their expected theoretical values from uniform sampling.

3) Partial Correlation Analyses

In Figure 2,

• The standard Pearson correlation would have certainly observed a significant correlation between variables \(A\) and \(C\).
• This spurious connection between \(A\) and \(C\) due to the confusing effect of \(B\) is clipped using partial correlation instead.

digraph G {
    layout=neato
    A -> B;
    B -> C;
}

3) Patient stratification

ii) Parameters’ estimation I: Methods

1. Define parameters’ configuration per endotype (observed concentration values, ODE definition model, …)

2. Run numerical optimisation algorithm to infer the values of free parameters to recover altered cytokine concentrations.

3. Mock \(TNF-\alpha\) blockade by fixing the total concentration to 0:

   1. In COPASI, you switch the species’ status from reactions to fixed.

   2. With CoRC, use setSpecies function with initial_concentration = 0 and type ="fixed".

absolute_parameter_per_cluster <- purrr::map(fit_experiments, function(cluster) {
  CoRC::clearParameterEstimationParameters(model = healthy_model)
  CoRC::clearExperiments(model = healthy_model)
  return(runParameterEstimation(
    parameters = fit_parameters,
    experiments = cluster,
    method = list(method = "LevenbergMarquardt",
                  iteration_limit = 500,
                  tolerance = 1e-6),
    update_model = FALSE,
    randomize_start_values = FALSE,
    create_parameter_sets = TRUE,
    calculate_statistics = FALSE,
    model = healthy_model
  )$parameters)
})

1

One parameters’ estimation per endotype configuration. The number of observations, aka concentrations’ species, must exceed the number of free parameters to infer ⚠️

2

clearParameterEstimationParameters and clearExperiments guarantee a fresh start for each parameters estimation (equivalent to pandas.deepcopy()).

3

runParameterEstimation is used for estimating the 10 free parameters allowed to vary.

4

All numerical optimisation approaches to reduce the cost function between the observed steady-state concentrations and the expected ones are listed in Optimization Methods. For reproducible analyses, you need to provide the method and numerical optimisation hyper-parameters, such as iteration_limit or tolerance.

3) Patient stratification

ii) Parameters’ estimation II: Main Reproducibility Hurdles

1. 🚫 No computer-readable datasets (only hard-coded pdf tables) ➡️ Manual typing is an error-proned process.

2. 🚫 Averaged differential expression inputs instead of raw, individual values, preventing from:

   1. Reproducing differential analyses reported in Table 7, Lo et al. (2016).

   2. Running a multivariate regression model, thus artificially underestimating the uncertainty of the parameters.

3. 🚫 No supplying of the parameters’ estimation approach:

   3. Optimisation method
   4. Hyper-parameters for evaluating numerical convergence
   5. Individual parameter weights

3) Patient stratification

iii) Parameters’ estimation III: Results

Simulated Fold Change variations of the cytokine and T-cell concentrations before and after TNF-\(\alpha\) blockade, per endotype. Blue bars represent the results of pre-treatment; red bars represent the results of post-treatment.

• In original Cluster 4 (reduced abundances of Th1 and Th2 cell subsets), all species’ concentrations are decreasing

• In reproduced Cluster 4, concentrations of IL-10, Tregs and TGF-\(\beta\) are increasing with respect to healthy state conditions, an effect magnified by blocking the production of IFN-\(\gamma\) ➡️ noting the positive feedback loop between these three entities, and the overall negative regulatory effect of Tregs on T cells concentrations, this result is biologically more meaningful.

3) Model Annotation

1. General annotations of the model: name, biological processes modelled, publication if any, curation authors, …

2. Annotate the models species (nodes of the ODE model) with established EBI-EMBL ontologies to enable cross-linking between databases:

   • UniProt for proteins (including cytokines)

   • Cell Line Ontology (CLO) for immune cells

   • Gene Ontology: GO terms for pathways.

3. (Optional): detail each ODE reaction, namely edges (for instance, provide bibliographic references mentioning interplay between two cell types).